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The effect of small scale forcing on large scale structures in /5-plane two- 
dimensional (2D) turbulence is studied using long-term direct numerical 
simulations (DNS). We find that nonlinear effects remain strong at all 
times and for all scales and establish an inverse energy cascade that ex- 
tends to the largest scales available in the system. The large scale flow 
develops strong spectral anisotropy: k~^/^ Kolmogorov scaling holds for 
almost all 4>, (p = aicta.n{ky/kx), except in the small vicinity of kx = 0, 
where Rhines's k~^ scaling prevails. Due to the k~^ scaling, the spectral 
evolution of /3-plane turbulence becomes extremely slow which, perhaps, 
explains why this scaling law has never before been observed in DNS. Sim- 
ulations with different values of /? indicate that the /3-effect diminishes at 
small scales where the flow is nearly isotropic. Thus, for simulations of 
/3-plane turbulence forced at small scales sufficiently removed from the 
scales where /3-effect is strong, large eddy simulation (LES) can be used. 
A subgrid scale (SGS) parameterization for such LES must account for 
the small scale forcing that is not explicitly resolved and correctly accom- 
modate two inviscid conservation laws, viz. energy and enstrophy. This 
requirement gives rise to a new anisotropic stabilized negative viscosity 
(SNV) SGS representation which is discussed in the context of LES of 
isotropic 2D turbulence. 



Preprint submitted to Elsevier Preprint 



5 February 2008 



1 Introduction 



Turbulent flows subjected to differential rotation develop spectral anisotropy. 
Understanding and modeling such flows present a major theoretical and exper- 
imental challenge, mainly due to their geophysical, astrophysical and plasma 
physics importance. The simplest two-dimensional system of this kind de- 
scribes the flow of a thin layer of homogeneous fluid on the surface of a ro- 
tating sphere. Such flows are governed by the barotropic vorticity equation in 
the /9-plane approximation in which the fluid moves in tangential planes 

Here, ( is the fluid vorticity, z/q is the molecular viscosity, and ^ is the forcing; 
X and y are directed eastward and northward, respectively. 

The constant (3 is the background vorticity gradient describing the latitudinal 
variation of the normal component of the Coriolis parameter, / = /o + (3y- 
Considering the energy transfer subrange typical of classical isotropic two- 
dimensional turbulence we assume that the forcing ^ is concentrated 
around some high wave number and is random, zero-mean, Gaussian and 
white noise in time with the correlation function 

(e(k, t) e(k', t')) = Stt r/ 5(P - kl) 5(k + k') 5it - t'), (2) 

where (. . .) denotes ensemble average. Thus defined, the forcing ^ supplies 
the system (1) with energy and enstrophy at rates / (|^(k, t)\'^)/k'^ dk/ (27r)^ = 
e = and / (|^(k, dk/{2iiY = rj, respectively. In the inviscid limit, the 
corresponding total energy and enstrophy vary in time as in the isotropic case, 
according to E(t) = E{0) + et and Q{t) = Q{0) + rjt; the /3-term does not 
enter these evolution laws explicitly. 

In the linear limit with no forcing, Eq. (1) describes propagation of planetary, 
or Rossby waves with the dispersion relation 

o'k = -Pkje, (3) 

while in the nonlinear limit without the /3-term, this equation represents clas- 
sical isotropic 2D turbulence. Thus, Eq. (|ID describes the interaction of 2D 
turbulence and Rossby waves and is important for understanding of planetary 
scale geophysical processes. 

The relative simplicity of Eq. ([l|) has placed it at the focal point of theoreti- 
cal geophysical fluid dynamics, and much has been learned from this model. 
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However, the large scale behavior of its solutions still remains controversial; 
their spectral evolution laws have not been well established, while the spec- 
tral anisotropy and the importance of nonlinearity have received insufficient 
attention in the existing literature. Here, using DNS, we address the prob- 
lems of scaling laws and spectral anisotropy of /3-plane turbulence governed 
by (1) and forced at small scales. In the next section, evolution and scaling 
laws of /9-plane turbulence will be discussed. On relatively small scales the 
/3-effect is small and the flow is nearly isotropic. On larger scales, the effect of 
the /3-term becomes progressively stronger, and spectral anisotropy develops. 
Along with that, spectral evolution slows down considerably, compared to the 
isotropic case. As a result, prohibitively long integration time is required for 
low wavenumber modes to become excited even for relatively moderate reso- 
lution. To address this problem, LES of /5-plane turbulence can be sought. In 
section 3, an eddy viscosity for such LES is calculated. This SGS representa- 
tion includes negative Laplacian term that accounts for the unresolved small 
scale forcing and inverse cascade. Thus, this SGS representation gives rise to 
a stabilized negative viscosity (SNV) SGS parameterization 0]. In section 4, 
some results of the application of this SNV representation to isotropic 2D flows 
are presented. Finally, in section 5, we summarize our results and conclusions. 



2 Evolution and Scaling Laws of /5-Plane Turbulence 

Preliminary asymptotic analysis of Eq. (1) reveals that as ^ oo, the ampli- 
tude of the /3-term decays as coscp [0 = aictaji{ky/kx)], so that for large 
k, the /3-effect is expected to be small and the system should behave mostly 
as isotropic 2D turbulence in the energy transfer subrange . The energy 

spectrum defined as E(k, t) = 2tt (|C(k, t)\^)/k, with / dk J^"" E(k, t)rf0/(27r) = 
E{t) being the total energy, is nearly isotropic, E(k,t) ^ E{k,t), and Kolmo- 
gorov-like, i.e. 



where Ck is the Kolmogorov constant. On the other hand, Rossby waves dom- 
inate as — >■ 0. Assuming that the energy spectrum in this limit is determined 
by P and k only, one finds, using dimensional considerations 



where Cr, by analogy to the Kolmogorov constant, can be referred to as the 
Rhines constant. For isotropic 2D turbulence, the estimate of the Kolmogorov 
constant Ck from DNS and analytical theories of turbulence is about 6. On 
the other hand, the Rhines spectrum (5) has not been observed in DNS so far 



EK{k) = CKe'/n-'/\ 



(4) 



ER{k)=CR(3^k-\ 



(5) 
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so that not only no estimate of Cr is available, but even the very existence of 
the spectrum (5) has been in doubt. 

Turbulence and waves have comparable effects when the eddy turnover time 
of isotropic 2D turbulence approaches the Rossby wave period, which occurs 
at a transitional wave number, kt{(f)), 



h{<j)) = kpcos'/'<j), kf,= {l3^/ef'\ (6) 

The contour (P) in wavevector space has been termed "the dumb-bell shape" 
by Vallis and Maltrud 0; earlier, a similar contour was called "lazy 8" by 
Holloway in the context of /?-plane flow predictability studies. According to 
0, as A; ^ kp, the inverse energy cascade slows down, and increasing spectral 
anisotropy facilitates preferential energy transfer into large scale structures 
with kx —>■ 0, or zonal jets. This conceptual picture was further developed in 
later publications. In particular, Vallis and Maltrud found that a.s k —>■ k^, 
the dumb-bell shape becomes a barrier for turbulence energy. However, since 
the dumb-bell shape excludes the axis k^ = 0, the inverse cascade can continue 
in the vicinity of kx = facilitating the funneling of the energy into coherent 
structures corresponding to zonal flows. Rhines 0, Vallis and Maltrud 0, 
and Holloway note that the spectral anisotropy induced by the /3-term can 
be associated with the mechanism of generation and maintenance of zonal 
flows. Zonal structures are robust features of the simulations by Rhines 0, 
Vallis and Maltrud [0] and Panetta [0, as well as of the large scale Earth 
and planetary circulations. On the other hand, for k < kp one would expect 
that nonlinear transfer becomes relatively ineffective, and the flow regimes 
inside the dumb-bell shape are more prone to linear dynamics Although 
this assumption underlines some theories, it has not been thoroughly tested in 
DNS so far. Generally, processes at k < kp were avoided in the previous DNS 
of /3-plane turbulence. Thus, not only could these processes not be properly 
resolved {kp usually did not exceed 20 or so), but the whole region k < kp 
was treated rather as an area of large scale energy dissipation due to linear 
(Ekman) drag or other factors. 

In the present DNS, the processes at scales of the order of k'^^ and larger 
are of primary interest. These DNS employ a fully-dealiased pseudospectral 
Fourier method to solve Eq. in a box {x, y} G [0, 27r] x [0, 2n] with doubly 
periodic boundary conditions and zero initial data. The spatial resolution is 
512^ and the 2/3 dealiasing rule is used. The time discretization is the same 
as in A hyperviscous term of the form z/5 (— 1)^"^^ A^C -|- (—1)'?"*"-^ A~'^C 
with p = 8 and g = 5 is introduced in (|T]) to ensure sharp energy removal 
from both small and large wave number modes. The forcing ^ is represented 
by = A^ak{t)/y/6i for k E [100,105] and is zero otherwise; here, 6t 

is the time step, and ak(t) is a discrete and uncorrelated in time, Gaussian, 
random number with unit variance. Parameter setting for DNS reported here 
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is: St = 1.0, ul = 20.0, z/5 = 1.0 x 10~^^, and = 0.1 which results in an 
energy injection rate of e ~ 2.7 x 10"^'^. Results of three DNS, with f3 = 0, 
0.053, and 0.3, corresponding to = 0, 56, and 158, respectively, are reported 
here and shown in Figs. 1-5. These cases will be referred to as Case 1,2, and 
3, respectively. All simulations were terminated when the energy reached the 
largest scales available in the computational box and condensation started 
affecting flow behavior at smaller scales |T3| , p!^ . 

In Figure 1, we plot the evolution of total energy and enstrophy normalized by 
the corresponding values, Eq and Qq, of the steady state, isotropic (/3 = 0) case. 
A striking feature is that increasing in (3 causes increasing in computational 
time necessary for the effects of condensation to appear. For example, for 
P=0.3, the effect of large scale energy saturation becomes significant at about 
t/Ttu ~ 50, where = 27r/ \/2Eq is the large scale eddy turnover time of 
the steady state in the isotropic case. For comparison, in the corresponding 
isotropic DNS, the large scale energy condensate forms after only about 4rf„. In 
fact, the dramatic increase in computational time necessary for the energy to 
reach the lowest active modes in the system is the major factor hampering DNS 
of /3-plane turbulence; the problem is exacerbated with increasing resolution. 
The reason why such a long integration time is required will be explained 
later. 

The results plotted in Figure 1 show that in agreement with the conservation 
law, E{t) exhibits linear in time growth. The analysis of mode- wise energy dis- 
tribution (not shown here) indicates that increasingly longer time is required 
for saturation of the lower modes. 

In Figures 2-4 we plot k-dependent energy spectra averaged over small sectors 
±7r/12 around = and = 7i/2. For Case 1, as expected, the spectrum is 
isotropic and Kolmogorov-like. When /? 7^ 0, spectral anisotropy develops. As 
Figs. 3 and 4 indicate, this anisotropy becomes more and more pronounced 
with increasing (3. Although, as mentioned earlier, spectral evolution slows 
down with /3 7^ 0, spectral anisotropy develops relatively fast. In Figure 3 
we show that for Case 2 with G [100, 105] and kj^ = 56, anisotropization 
begins at t ~ r^^, while for Case 3 with kj^ = 158, anisotropization starts 
immediately. Long term integration indicates that in the direction = 0, or 
ky = 0, the Kolmogorov scaling (4) develops, while for the directions = 
ib7r/2, a scaling similar to that of Rhines (5), k~^, prevails. Compared to Fig. 
3, Fig. 4 presents results of DNS with larger P, or larger k^, leading to better 
resolution of the region k < kp and a more pronounced k~^ spectrum. In 
Figure 5, we compare the compensated spectrum E{k)(3~'^k^ along = ±7r/2 
at the end of integration for Cases 2 and 3, i.e., for t/rtu = 16.25 and 30.12, 
respectively. An important conclusion that can be drawn from this comparison 
is that along = ±7r/2, the energy spectrum scales with Combined with 
dimensional considerations, this result reaffirms the validity of the scaling (5) 
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along = ±7r/2. More detailed analysis of the anisotropic spectra reveals 
that the Kolmogorov scaling holds for almost all directions except 

for the small vicinity of = ±7r/2. It is important to emphasize that the 
inverse energy transfer does not cease for any direction and any k < kp and 
eventually extends to the largest scales available in the system. Rather than 
being a barrier, or even a soft barrier for the inverse energy cascade, kp only 
appears to be a threshold of spectral anisotropy. After the energy front reaches 
kp, some readjustment of Ck is observed; its value decreases from 6 to about 3. 
This points to the active energy exchange between the Kolmogorov and Rhines 
scaling regions facilitated by nonlinear interactions and to intensification of 
the Rhines's flow regime at the expense of its Kolmogorov counterpart. 



Generally, since the k~^ spectrum is much steeper than the k~^^^ spectrum, 
most of energy resides in sectors closely adjacent to = ±vr/2, so that the 
angular-averaged energy spectrum also obeys Rhines's k~^ scaling (5). There- 
fore, the present DNS confirms Rhines's arguments in 0. Indeed, this is the 
first time that the k~^ spectrum is observed in DNS. Our results also show 
why meaningful DNS of /5-plane turbulence requires much longer integration 
times than the corresponding isotropic DNS. Indeed, from the total energy 
evolution law for both isotropic and /3-plane turbulence, E{t) oc et, one es- 
timates that the time required for the energy front to reach wave number k 
is t oc E{t)/e. For the same e, the ratio of these characteristic time scales, 
tn/tK, for Rhines and Kolmogorov spectra can be estimated from (4) and (5) 
by angular averaging and integrating corresponding energy spectra between 
oo and k giving tji/tx oc {kp/kY^^^. For Cases 2 and 3 and k = 15, one finds 
that tn/tx is of the order of 80 and 2500, respectively. 



Here, further discussion of the physical meaning of kp given by (6) would be 
relevant. Other definitions of fc^ have been used in the literature interchang- 
ingly with (6) (see, for instance, 0). One of those definitions is due to Rhines 
H, = {(3/2UY'^, and another one is due to Holloway and Hendershott |15| 



kp = (3/C- Here, f/ is a characteristic velocity which can be identified with 
the square root of the total energy of the system, while ^ is a characteristic 
vorticity that can be identified with the square root of the total enstrophy of 
the system. For the k~^ spectrum, it can be shown that, within numerical co- 
efficients, both Rhines's and Holloway and Hendershott's definitions provide 
kj3 = k, where k is the smallest wavenumber attained by the energy front. 
Furthermore, taking into account the evolution law of /3-plane turbulence, one 
can show that thus defined kp is time dependent, kp ~ Therefore, defi- 

nitions (6) and those suggested by Rhines ||^ and Holloway and Hendershott 
||T5| characterize different processes in /3-plane turbulence. While the former 
provides the stationary threshold separating the regions of nearly isotropic and 
strongly anisotropic, Rossby wave dominated turbulence, the latter pertain to 
the largest scales available in the system at any given time. 
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To better understand the anisotropy of the energy (or enstrophy) spectrum 
of /5-plane turbulence, let us consider the enstrophy transfer function 7n(k, t) 
derived from the enstrophy equation 



!J(k,i) =rn(k,i). (7) 



Here, fi(k, t) is vorticity correlation function, and 7Q(k, t) is given by 

T^(k,t) = 2vrA: J ^(C(p, t)C(q, t)C(-k, t))^. (8) 

p+q=k 

Let us introduce a cutoff wavenumber kc separating explicit {k < kc) and 
implicit, or SGS {k > kc) modes and consider enstrophy exchange between all 



implicit modes and a given explicit mode k in the limit t ^ oo. Following |16 
this exchange is denoted by '/^(klfcc) and is obtained from (8) by extending 
integration only over all such triangles (k, p, q) that p and/or q are greater 
than kc- 

In Figure 6, we plot the anisotropic energy transfer 7£;(k|/cc) derived from the 
corresponding enstrophy transfer 7Q(k|A;c) for Case 3 with kc = 50. Obviously, 
TE(k\kc) describes energy transfer from all SGS modes with k > kc to an 
explicit mode with the wavenumber h, k < kc- As in |T2|, this energy transfer 
function was calculated here directly from DNS results. Consistently with the 
isotropic case |]T2[, TE(k\kc) develops a cusp at A; — * kc. As A; — > 0, strong 
anisotropy prevails, and most of the energy is funneled into the sectors adja- 
cent to = ±7r/2. On the one hand, the modes with small fc^ carry most of the 
energy; on the other, they allow for weak x-dependency of the flow field, thus 
maintaining nontrivial nonlinearity which in turn sustains anisotropic transfer. 
Thus, the modes with small k^ correspond to nearly one-dimensional, zonal 
structures, or jets. Such structures are clearly seen in an instantaneous snap- 
shot of the vorticity field in physical space at t = 39.6Ttu for Case 3 as plotted 
in Fig. 7; these jets have also been a robust feature of earlier simulations |l7|,p!0[|. 
The observation that can be made based upon the present DNS is that the 
number of jets is very close to the minimal excited wavenumber in the energy 
spectrum; there is no obvious linkage between the number of jets and fc/j. 

The zonally- averaged velocity component U {y, t) is plotted in Fig. 8 for Case 3. 
Initially, the jets were nearly symmetric with respect to the reflection y — * —y. 
Bjt = 39.6Ttu, the zonal jets develop strong asymmetry; as also observed in ||^, 
eastward jets are sharp and narrow while westward jets are smooth and wide. 
It is important to know whether or not these jets are stable with respect to 
perturbations with nonzero k^. Although the velocity field is time-dependent, 
the meridional motions are rather slow compared to zonal flows. In this case. 
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the Rayleigh-Kuo inviscid stability criterion generalized for /3 7^ requires 
that the profile U{y,t) — (3 y"^ /2 has no inflection points for the flow to be 
linearly stable. Examination of the second derivative Uyy{y,t) shown in Fig. 
8 demonstrates that the Rayleigh-Kuo criterion Uyy{y,t) — /5 7^ does hold, 
which also agrees with the results in Furthermore, as time advances, the 
profile of Uyy{y,t) develops large negative peaks thus ensuring linear stability 
at later times. As long as the effect of the large scale drag remains small, 
excitation of progressively smaller wavenumber modes leads to a diminishing 
number of jets, in agreement with other simulations It is interesting to 
note that in models of barotropic flows on a /3-plane with externally imposed 
mean flow U, the difference (3e = 13 — Uyy may decrease to zero or even become 
negative. The sign of j3e is related to interaction between the critical layer 
and Rossby waves and may point to absorption {(3^ > 0), perfect reflection 
{(3e = 0), or over- reflection (/?e < 0) of Rossby waves by the critical layer 



(for more details, see |[T7|-pT|). In the present DNS, as well as in the previous 
studies by Vallis and Maltrud 0, the mean flow U is generated and sustained 
by small scale random forcing, while (3e is always positive. 



In summary, it is useful to highlight the peculiarities caused by the nonlinear 
terms in barotropic vorticity equation on a /3-plane. Although the /3-term does 
not enter the energy and enstrophy equations explicitly, it has a profound effect 
on the energy spectrum and spectral transfer. This effect is solely due to the 
nonlinearity of Eq. (1) which facilitates interaction between Rossby waves 
and vorticity modes. As a result, the energy spectrum and transfer develop 
strong anisotropy for k < kp while the inverse energy cascade extends to 
ever smaller k. Furthermore, inside the dumb-bell shape @, where, by simple 
scaling considerations, the /9-effect is expected to prevail, the energy spectrum 
is determined not by j3 but by the presumably irrelevant parameter e. On the 
other hand, in the small sectors around = ±7r/2 outside the dumb-bell 
shape naive considerations show that the /3-term is vanishing [j3dxV~'^C —>■ 
il3kxk~'^C, ^ as fca; — > in (1)] and the mechanism of anisotropic inverse 
transfer funnels energy into zonal jets. In the present DNS we find that the 
zonal jets do indeed form, but the effect of the /3-term does not diminish by 
any means: the spectrum of energy in the regions — > ±7r/2 is determined by 



The possibility of energy transfer from Rossby waves to zonal flows has long 
been discussed in the literature and by no means is a trivial matter. Straight- 
forward generalization of the Fjortoft theorem |^2[ shows that a resonant triad 
of wavevectors one of which is aligned with the axis y does not allow for direct 
energy transfer into zonal flows. However, the vector with fc^. = facilitates 
energy exchange between the other two vectors - members of the triad. This 
fact was proved in in much more general way; it was shown that the cou- 
pling coefficient is zero for any such triad that kx = for one of the vectors. 
Later, Newell EM showed that Rossby waves can transfer energy into zonal 
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flows through a sideband resonance mechanism or through a quartet reso- 
nance. Reinforcing this result, the present DNS indicate that not only energy 
transfer into zonal flows due to nonlinear interaction is possible, but in fact it 
is the dominant mechanism of /3-plane turbulence on large scales. 



3 Parameterization of Eddy Viscosity for LES of /?-Plane Turbu- 
lence 



Since DNS of /9-plane turbulence in the energy transfer subrange demands ex- 
cessive computational resources even for relatively moderate resolution, LES 
could be a viable alternative. For such LES, however, a SGS representation 
is crucial since not only should it accommodate inverse energy transfer but 
also it should account for the small scale energy forcing excluded in the ex- 
plicit modes. Following [Q, an SGS representation for high wavenumber forced 
/5-plane turbulence can be given by an anisotropic two-parametric viscosity, 
z/(k| /cc), defined by 

n\u\ ^n(k|A;c) 



This two-parametric viscosity was calculated analytically using the renormal- 
ization group (RG) theory of turbulence |]2S|; it is shown in Fig. 9 for the set 
of parameters corresponding to Case 2 of DNS. The behavior of the theoret- 
ically derived z/(k|fcc) is quite similar to that obtained from DNS as shown 
in Fig. 10 also for Case 2. For large k, the /3-effect is weak, and i/(k|fcc) be- 
haves similarly to the isotropic case |25|JT^ ; there is a sharp positive cusp and 
then i/(k|/cc) becomes negative. As — 0, the effect of the /3-term becomes 
stronger; z/(k|/cc) remains negative in the vicinity of = ±7r/2 but increases 
to zero in other directions. The negativity of i/(k|fcc) along = ±7r/2 is con- 
sistent with the strong energy flux into zonal jets as discussed in the previous 
section. It is important to mention that when the DNS integration was ter- 
minated, the low wavenumber end of the spectrum was not yet fully energy 
saturated. This fact explains the sharp increase in DNS-derived ^(klkc) at 
small k evident in Fig. 10. The sign-changing structure of i/(k| /cc) is indicative 
of a physical space SGS representation which combines a negative Laplacian 
and a positive (dissipative) hyperviscosity thus presenting an anisotropic ver- 
sion of the Kuramoto-Sivashinsky equation or SNV SGS representation 0]. 
The main ideas of the SNV formulation are discussed in the next section as 
applied to LES of more simple isotropic 2D turbulence. The isotropic SNV 
SGS representation can be useful for LES of /3-plane turbulence with fc^ < kc 
since deviation from isotropy is expected only for k < kp where, due to the 
factor of k"^ on the right hand side of (10) below, the SGS contribution is small. 
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4 SNV SGS Representation for Isotropic 2D Turbulence 



Specific difficulties of LES of 2D turbulence in the energy transfer subrange 
are that the SGS representation must simultaneously be consistent with the 
conservation laws for both energy and enstrophy (or potential enstrophy) 
In addition, in LES of the inverse energy cascade, the wavenumber modes 
of the energy source lie in the subgrid scales so that the SGS representation 
must assume the function of the forcing. The SGS schemes being used to date 
in both 3D and quasi-2D flows usually employ Laplacian or higher order hy- 
perviscosities that ensure SGS energy dissipation but not forcing. While such 
SGS are consistent with the dynamics of 3D turbulence where energy is di- 
rectly cascaded from large to small scales where it dissipates, they cannot be 
expected to work well in quasi-2D flows where energy cascades from small 
to large scales while the enstrophy is transferred in the opposite direction. A 
possible resolution to this problem can involve the replacement of the SGS 
forcing by force located in the explicitly resolved region near kc- However, 
this solution is not only quite cumbersome but also significantly distorts the 
explicit scales near kc- In addition, this approach is difficult for implementa- 
tion in physical space, particularly for bounded systems and/or systems with 
spatially nonuniform energy sources. 

An alternative approach is offered by the SNV representation developed in 
1^ for isotropic 2D turbulence. In this approach, Eq. (1) (with f3 = 0) that 
includes small scale forcing ^ is replaced by LES equation 

^+ / ^<(pm-p)f^ = -m.)k'm, do) 

\p\,\k~p\<k, 
< k < kc- 



Here, eddy viscosity is identified with the two-parametric viscosity z/(A;|A;c) 
defined by (9) (where, due to isotropy of the fiow, the angular dependence is 
omitted). In it was found that this eddy viscosity is given by 

^ik\kc) = 77^N{k/k,), (11) 



where Q{t) is the total enstrophy of the system and N{k/kc) = ^{klkc) /\i'{0\kc)\ 
is plotted in Fig. 11. This formulation is designed to provide constant energy 
transfer with the rate e. Consistent with the eddy viscosity approach, the SGS 
representation (11) is a function of the fiow. In a statistical steady state, Eq. 
(11) coincides with the two-parametric viscosity derived from the RG theory 



10 



In a series of LES that employ the formulation (11) together with a specially 
designed large scale drag that effectively removes energy that reaches the 
largest computational scales, it was found that simulations can be carried out 
nearly indefinitely (being terminated after about 200 turnover times), while a 
very robust Kolmogorov k^^^^ energy spectrum was established. 

The formulation (11) can be simplified to make it more easily adaptable for 
simulations of quasi-2D flows in physical space. First, a two-term power series 
expansion in powers of k/kc that satisfies all necessary conservation laws can 
be derived from (11) 0], 

v{k\kc) 



25 



l8^](t) 



-f-V 

5 \kj 



(12) 



Here, the negative term in the square brackets accounts for the effect of the 
small scale forcing and inverse energy cascade. The second term describes dis- 
sipative processes and is equivalent to a biharmonic viscosity in physical space. 
Then, (12) can be further simplified if Vt in its dissipative term is replaced by 
the value obtained from the Kolmogorov spectrum (4). Finally, the resulting 
SGS representation for 2D turbulence in the energy transfer subrange takes 
the form 

^{k\K) = + 0.511eV3fc;io/3^2_ (^3) 



In Figs 12 and 13, we plot the results of LES of 2D turbulence with the SGS 
representation (13). One can see that the total energy and enstrophy oscillate 
slightly around their steady state values while a robust Kolmogorov energy 
spectrum is established. 

The corresponding physical space representation of the SGS operator for quasi- 
2D flows in the energy subrange is obtained using the inverse Fourier transform 
of (13), 




where 



^4 = 0.511e^/3(A/27^)l°/^ (16) 
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and where 17 (x) denotes the enstrophy averaged over an area adjacent to 
the grid cell, A is the grid resolution, the Laplacian term in (14) is written in 
the conservative form. Since (14) includes two terms, a negative Laplacian and 
positive (in the sense of dissipation) biharmonic, it is referred to as a stabilized 
negative viscosity (SNV) formulation. Generally, SNV equations are far more 
complicated than equations of the Kuramoto-Sivashinsky type p6| , ^ because 
in the former, the coefficients are not constant but, as in the classical eddy 
viscosity approach, are functions of the flow. 



5 Conclusions 



In this paper, we describe results of direct numerical simulations of /?-plane 
turbulence forced at small scales. As for isotropic 2D turbulence, nonlinear 
interactions play a crucial role in the evolution and dynamics of /?-plane tur- 
bulence on all scales and at all times. The effect of the small scale forcing 
propagates to ever larger scales due to the inverse energy cascade. On rel- 
atively small scales, k > k[^, the flow is nearly isotropic and behaves like 
classical 2D turbulence. For modes with eddy turnover times comparable with 
the Rossby wave period, the turbulence energy spectrum begins to develop 
anisotropy that increases for — > 0. It is found that for almost all (j) except 
for the small vicinity of = ±7r/2, the spectrum of turbulence is Kolmogorov- 
like, E(k) ~ k~^^^. Around (p = ±7r/2, or \k.j.\ — > 0, the two-parametric en- 
ergy transfer function TE(k\kc) attains a sharp maximum, which indicates 
that there is a strongly anisotropic energy transfer into quasi-one-dimensional 
(zonal) structures. These zonal flows are generated and sustained by nonlinear 
processes, particularly, anisotropic inverse energy cascade. For — ±7r/2, the 
energy spectrum becomes very steep, E{k) ~ k~^, as was first suggested by 
Rhines 0. The steepening of the spectrum slows down fiow evolution, to the 
extent that DNS even with relatively moderate resolution becomes computa- 
tionally prohibitive. This difficulty can be resolved through LES of /5-plane 
turbulence. For such LES, an SGS representation has been designed based 
upon a two-parametric viscosity. We have applied this approach for an exam- 
ple of isotropic 2D turbulence for which case the SGS representation is given 
by the SNV parameterization. 
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Fig. 1. Evolution of total energy E{t) (top) and enstrophy n{t) (bottom) for Cases 
1, 2 and 3 (solid, dotted, and dashed lines, respectively). 
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Fig. 2. Energy spectra for = (dotted line) and (p = 7i"/2 (dashed line) averaged in 
time and over a small surrounding sector zb7r/12 for t/rtu = 0.63, 1.19, and 8.32 for 
P = (Case 1) (due to the central symmetry, only these two directions are shown). 
Long dashed lines depict slopes of —5 and —5/3. 
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Fig. 3. Same as in Fig. 2 but for t/rt„ = 0.63, 1.19, and 16.25 and /3=0.053 (Case 
2). 
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Fig. 4. Same as in Fig. 2 but for tjTtu = 0.63, 1.19, and 30.12 and /3=0.3 (Case 3). 




Fig. 5. Compensated energy spectrum for Cases 2 (dashed line) and 3 (solid line) 
for t/rtu = 16.25 and 30.12, respectively. 
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Fig. 6. Spectral energy transfer, TE(k|A;c), for kc = 50. 

Fig. 7. Instantaneous vorticity field C,{'x.,t) at t/Ttu = 39.6. [Too hig for this archive, 
available upon request.) 
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Fig. 8. Zonally-averaged velocity profile U{y, t) (left column) and its second deriva- 
tive Uyy{y,t) (right column) for Case 3 at (from top to bottom) t/Ttu= 21.8, 39.6, 
and 59.4. 



21 




Fig. 9. RG-derived two-parametric viscosity normalized by the absolute value of its 
minimum for /3-plane turbulence. Parameter setting corresponds to Case 2 DNS. 
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Fig. 10. DNS-based two-parametric viscosity normalized by the absolute value of 
its minimum for Case 2 of /3-plane turbulence. 
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Fig. 12. Evolution of the total energy E{t) (a) and enstrophy ^(t) (b) in LES of 
isotropic 2D turbulence. In (a) we also plot the evolution of E{t) with the energy 
of the 4th, 5th, 6th and 7th modes removed. 
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Fig. 13. Time averaged energy spectrum for LES of isotropic 2D turbulence. The 
solid line depicts the Kolmogorov —5/3 slope. 
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